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ABSTRACT 

Context. The X-ray spectra observed in the persistent emission of magnetars are evidence for the existence of a magnetosphere. The 
high-energy part of the spectra is explained by resonant cyclotron upscattering of soft thermal photons in a twisted magnetosphere, 
which has motivated an increasing number of efforts to improve and generalize existing magnetosphere models. 
Aims. We want to build more general configurations of twisted, force-free magnetospheres as a first step to understanding the role 
played by the magnetic field geometry in the observed spectra. 

Methods. First we reviewed and extended previous analytical works to assess the viability and limitations of semi- analytical ap- 
proaches. Second, we built a numerical code able to relax an initial configuration of a nonrotating magnetosphere to a force-free 
geometry, provided any arbitrary form of the magnetic field at the star surface. The numerical code is based on a finite-diff'erence 
time-domain, divergence-free, and conservative scheme, based of the magneto-frictional method used in other scenarios. 
Results. We obtain new numerical configurations of twisted magnetospheres, with distributions of twist and currents that diff'er from 
previous analytical solutions. The range of global twist of the new family of solutions is similar to the existing semi- analytical models 
(up to some radians), but the achieved geometry may be quite diff'erent. 

Conclusions. The geometry of twisted, force-free magnetospheres shows a wider variety of possibilities than previously considered. 
This has implications for the observed spectra and opens the possibility of implementing alternative models in simulations of radiative 
transfer aiming at providing spectra to be compared with observations. 

Key words, magnetohydrodynamics - stars: neutron - stars: magnetic fields - scattering 



1. Introduction 

Anomalous X-ray pulsars (AXPs) and soft gamma-ray re- 
peaters (SGRs) are a class of neutron stars (NSs) character- 
ized by high X-ray quiescent luminosities, short X-ray bursts, 
and giant flares (for SGRs). They are believed to be magne- 
tars (iDuncan & Thompson. 1992, Thompson & Duncan, 1996): 
young (P/2P ~ 10^-10^ yr) NSs with very strong magnetic 
fields ~ 10^^ - 10^^ G. Magnetars are supposedly born with 
a very short spin period (few milliseconds), which causes the 
formation of an intense inner magnetic field via dynamo pro- 
cesses and a subsequent rapid loss of rotational energy by means 
of magnetic braking. As more and better data become available 
in the last decadeQ the separation between the two traditional 
classes (SGRs and AXPs) has become thinner, and they are no 
longer considered different classes. All SGRs and most of AXPs 
have shown sporadic bursts in the X and y bands, and the recent 
discovery of transient magnetars has made the situation even less 
clear. Currently, our interpretation of the term magneta rs is be ing 
reconsidered (see Mereghetti (20081). iRea & Esposito (l201ll) for 
observational reviews). In general, magnetars are characterized 
by long periods P ~ 2 - 12 s and high X-ray persistent luminosi- 
ties Lx ~ 2 X 10^^ - 2 X 10^^ erg s"\ which point to the magnetic 
field decay as the main source of energy, unlike radio pulsars, 
fed by the loss of rotational energy. One key ingredient in those 
models that try to explain the variety of magnetar phenomenol- 
ogy and their spectra is the presence of toroidal fields, of the 
same order of magnitude or even larger than the poloidal com- 
ponent, both in the interior and in the external magnetosphere. 



^ The updated catalog of these sources can be found at 
http : //www . physics . mcgill . ca/~pulsar /magnetar /main . html 



The object of this paper is the modeling of the magneto- 
sphere of magnetars, which is crucial because it is one of the 
main ingredients in the formation of spectra. The difl&culty in 
building fully consistent solutions of NS magnetospheres is such 
that, after more than 40 years, the aligned rotator model pro- 
posed by iiGoldreich & JulianI 11969) still remains very popular 
for explaining the basic features of a rotating magnetosphere. We 
must point out that the aligned rotator model, with all its variants, 
is surely a naive description of pulsars that has been strongly 
criticized (e.g. Michel, 1980), among other reasons, simply be- 
cause it cannot explain why pulsars pulsate. A further refine- 
ment consists of the ideal MHD approximation for force-free 
magnetospheres, which leads to the so-called pulsar equation. 
This is simply the balance of the electromagnetic forces around 
the rotating neutron star under the assumption of axial symme- 
try and it neglects the pressure, inertial, and gravitational terms. 
The rotationally induced electric field is required to be perpen- 
dicular to the magnetic field. The latter can in general have an 
azimuthal component, that twists the poloidal field lines. This 
toroidal component has to be sustained by currents, whose shape 
is the source term of the pulsar equation: it is the fundamental 
parameter to be chosen in order to find the corresponding geo- 
metric configuration of the magnetic field. 

The problem is analytically solvable onl y for a f ew sina - 
ple, ad-hoc choices of the current, such as in iMichell (Il973bh . 
but these solutions generally present unphysical features in 
some border regions. A n alternative numerical approach by 
"Contopoulos et al.' ("1999) is useful to construct consistent force- 
free models, with a smooth matching across the light cylinder. In 
this model the corotating region only hosts an untwisted dipolar 
magnetic field, and the focus is on the bundle of open twisted 
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lines that cross the Hght c ylinder and the re sulting output power. 
The model was refined by GruzinoyI (l2QQ5l) and extended to the 
time-dependent, misali gned case by means of 3D simulations b y 
[Spitkovskvl(l2QQ6h and lKalapotharakos & Contopoulosl(r2QQ9h . 

Nowadays, X-ray instruments allow us to obtain more and 
more information from the spectra of magnetars. An especially 
interesting feature in the persistent emission of all of the magne- 
tar candidates is that their spectra can be well fitted with a ther- 
mal component (0.4 - 0.7 keV) plus a hard nonther mal tail, de- 
scribe d by a power law with photon index ~ 3 - 4 (iMereghettiL 
|2008l). This hard component is commonly attributed to resonant 
Compton scattering, which becomes important in the presence 
of a strong current (i.e. high charge density). For this reason, 
much attention has been given to the modeling of currents flow- 
ing in a twisted magnetosphere and obtaining estimates of the 
charge carrier density and the corresponding resonant optical 
depth. Lyutikov & Gavriil (2006) propose a semi- analytical ID 
model able to describe the basic properties of radiation trans- 
fer across a twisted magnetosphere. More accurate 3D Monte 
Carlo simulations have been performed for t he non relativistic 
(iFernand ez & Thompson. "2007", 'Nobili et al.', 2008a) and rela- 
tivistic cases (Nobili et al., 2008b). Fernandez & Davis (2011) 
discuss the implications of the outgoing polarization. The im- 
plementation in XSpec and the systemati c fits to observat ional 
data have been su ccessfully performed by lRea et"aD (l2008l) and 
IZane et al ] (12009) for the ID and 3D cases, respectively. All 
these very remarkable results rely on the sam e family of mod- 
els: th e self-similar twisted dipole propose d by [Thompson et al.1 
(l2002h . extended to higher multipoles by iPavan et al.l (12009'). 
which represents a semi-analytical, easy-to-implement solution 
to the pulsar equation. However, this model lacks generality, as 
it relies on a very particular choice of the form of the current. 

We note that force-free magnetosphere models rely on a 
purely macrophysical approach but a fully coherent microphys- 
ical description is still lacking: how are the currents and electric 
fields in a given model generated and sustained? Which parti- 
cles are present (electrons, ions, pairs)? Where do they come 
from? Which velocity distribution can be assigned to each kind 
of particle? In the magnetar framework, there are some excel- 
lent attempts to give a microphysical description of coronae. 
iBd oborodov & Thompson ( 2007) deal with the problem con- 
sidering a deviation from force-free equilibrium, which has to 
be introduced to explain the production of pairs, but many ques- 
tions remain open. Given the complexity of the problem, one 
has to decide whether to focus on detailed microphysical pro- 
cesses of an approximate macrophysical solution or to improve 
the large-scale electrodynamics solutions by gradually removing 
some of the simplifying hypotheses. 

In this article, we present our eff'orts to obtain some more 
general configurations of a twisted magnetosphere and to esti- 
mate the dependence of the spectrum on the geometric config- 
uration. In Sect. [21 we set the formalism needed to describe the 
problem and present some analytical and semi-analytical gener- 
alizations of previous works. In Sect. [3] we describe the numer- 
ical code used to build force-free configurations and the battery 
of tests performed to verify the code. We pay special attention 
to the eff'ect of key input parameters and the convergence of the 
method. In Sect. |4] we present new results for magnetosphere 
models, and explore the sensitivity of the relevant output (cur- 
rent, charge distribution, and resonant optical depth) to diff'erent 
input parameters. In Sect. [5] we summarize our findings and out- 
line future improvements. 



KP) 




Fig. 1. Relation between magnetic flux function F, enclosed cur- 
rent function /(F), and toroidal field. 

2. Analytical and semianalytical solutions 

2.1. The pulsar equation 

We consider an NS with a magnetic field B and angular velocity 
n. Assuming that the plasma is a perfect conductor in vicinity 
of the star, rotation induces an electric field given by 



E 



O X r 



xB. 



(1) 



where c is the velocity of light. Neglecting centrifugal (i.e. in- 
ertial), collisional, and gravitational terms compared to electro- 
magnetic forces, the force equilibrium equation is expressed by 



PeE -\--J xB = , 

c 



(2) 



where pe is the charge density and / the current density. 
Assuming axial symmetry (aligned rotator) and working in 
spherical coordinates (r, ^, 0), we introduce the magnetic flux 
function F(r, 0), in terms of which the poloidal field is 



^ VT(r,e)X(f> 
t>p = , 

r sin ^ 



(3) 



where is the unit azimuthal vector. With the choice of gauge 
F = on the axis, the function T(r,6) is related to the 
0-component of the potential vector by 



T(r,6)=A^(r,6)r sin 6, 



(4) 



and it is constant along a field line (Bp • VF = 0). Thus its value 
labels the axisymmetric surface Sr given by the azimuthal ro- 
tation of one field line. The magnetic flux flowing within the 
surface Sr, with value Or(^) = 27rF, is conserved by definition. 
A toroidal component (i.e. 0-component in axial symmetry) of 
the magnetic field is present if sustained by poloidal currents. 
The enclosed current / flowing within Sr (see blue arrows in 
Fig. [Hand note that the force-free condition implies that currents 
flow parallel to the surfaces Sr) has to be a function of only F, 
with the only requirement being that 7(0) = 0. Thus we have the 
freedom of choosing the enclosed current function /(F), and by 
Ampere's law, the 0-component of the magnetic field is related 
to /by 



cr sin 6 



/(F). 



(5) 



For a given magnetosphere model, it is useful to define the 
twist of a field line as the azimuthal displacement between its 
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surface footprints (for a dipolar-like configuration, one foot- 
print in each hemisphere). The integral of the local displacement 
S(f)(r,6) = B(pS6/Besm6 along the magnetic field line, /r, be- 
tween its surface footprints at magnetic colatitudes ^1,2 (r) is the 
line twist 



Be(r(0, T),0) sinO 



do. 



(6) 



where the dependence r(0, T) can be found by solving the field 
line equations. The function Acpt^i^) is used to characterize a 
magnetosphere. 

We then return to the equilibrium equation In axial sym- 
metry, only poloidal components of the electric field are al- 
lowed. Considering the poloidal component of Eq. ^ in cylin- 
drical coordinates (p, 0, z) and denoting the radius of the light 
cylinder by r/ = c/Q, we obtain the so- called p ulsar equation 
(IScharlemann & Wagoned ITOTllMichelL 1 1973bh : 



(7) 



where we indicate the partial derivatives of T with a subscript. 
The terms {plrif arise from the Coulomb force peE. The pulsar 
equation is actually the Grad-Shafranov equation for force-free 
fields. In general, Eq. d?]) has to be solved numerically. Even 
for trivial choices of untwisted configurations (/(F) = 0), the 
corotation velocity of the charged particles distorts the vacuum 
solutions, r esulting in op en field lines near the light cylinder 
p ~ 0{ri) (iMichei Il973 a). The only fully analytical solution 
including rotation and a smooth matching acro ss the light cylin - 
der is the split twisted monopole presented by'Michell (Il973bl) , 
in which the field lines are radial and the magnetic field changes 
sign between the northern and southern hemispheres. A current 
sheet at the equator is needed to ensure the divergence-free con- 
dition. This is a wind-like solution inappropriate to describing 
the region with closed field lines. Later, the numerical solution 
by Contopoulos et al.' (| 1999 |) for the first time provided a com- 
plete description of a magnetosphere, whose closed lines are un- 
twisted, while the wind zone has a configuration that is very sim- 
ilar to the split monopole. 

In the case of a nonrotating star (magnetars are slow rotators 
and this may be a reasonable approximation for some problems), 
no electric field is induced and the equilibrium equation is sim- 
ply JxB = 0: the force-free condition only requires that currents 
flow parallel to magnetic field lines. In this limit, we can refor- 
mulate the whole problem with the simple equation 



VxB = a(r)B , 



(8) 



where the function a(T) is related to the enclosed current func- 
tion introduced in Eq. ([5]) by 



2d/(F) 
c dF 



(9) 



The pulsar equation in spherical coordinates is now expressed as 



d^Y cos 6/ 1 ^F 1 d^Y 
^ ~ sin^^^ 



Next we expand the magnetic flux function in Legendre polyno- 
mials, Piifi): 



Y(r,ji) = Yoy-ai(r)(l-ji^) 

r ^* 



where ji = cos 0, ai(r) is the dimensionless radial function, 
the radius of the star, and Fq the magnetic flux normalization. 
Denoting the strength at pole by ^o, we hereafter choose 



ro = 



Borl 



(12) 



This leads to the following expression for the poloidal magnetic 
field components: 

Br = ^-V /(/+l)W^/(r), 



Bor* 
' 2 r 



V dPiilS) d(rai(r)) 
Z 



(13) 



d/j, dr 

We can obtain the governing diff'erential equation from Eq. ([TOb 



a{Y) a{Y)&Y . 



(14) 



Finally, by using the orthogonality relations of Legendre poly- 
nomials, we have 

Kl+\) \ d\rai{r)) ai{r) 
Bori.— — - r-^ /(/+ 1) 



2/ -hi I dr^ 

-1 



— — \a(Y) ai 
J-i djd [ J 



(r)dF 



dA/. 



(15) 



The righthand side depends on the functional form of a(Y) and 
is responsible for the coupling between difl'erent multipoles. For 
constant a the equation has analytical solutions (see Sect. 12.2b , 
while the problem is more difl&cult for other choices. The rest of 
this section is devoted to reviewing some possible choices that 
make the problem (semi)analytically solvable in the nonrotating 
case, which applies to magnetars. 

2.2. Constant a: Bessel solutions 

The choice of constant a = klr^, where ^ is a dimensionless pa- 
rameter, leads to decoupled equations for each single multipole, 
so that for each /, we must independently solve the equation 



2d^ai(r) dai(r) 
r h 2r — h 



dr^ 



dr 



k—\ -/(/+1) 



ai{r) = 0. 



(16) 



The analytic solutions of this equation are the spherical Bessel 
functions of the first and second kinds: 



Ji{x) = i-xY 



ld_ 

X dx 



smx 



FKx) = (-x)'-M-|- 

\xdx 



cosx 



(17) 



(18) 



where x = kr/r^. The constant ratio /k gives the typical scale- 
length of the magnetic field. The physical solutions (from which 
the vacuum solution Bp oc r"^^^^^ can be recovered in the limit 
^ ^ 0) are represented by the functions of the second kind, ex- 
plicitly: 



ai(r) = ciYi(x) , 
(11) Br=Y-Tj^(i^^)Pi(j^)^iYi(x)- 
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dPiiji) d(xYi(x)) 
^ 



(19) 



where c/ is the weight of the /-multipole. 

This choice of con stant a has been followed in several st ud- 
ies of the solar corona (Chiu & Hilton, 1977", "Seehafer, 1978) or 
applied to the region with o pen field lines of the pulsar magneto- 
sphere (IScharlemann & Wagoneiill973i) . The Yi{x) functions are 
oscillatory at large distances, and the magnetic field components 
change sign as r varies (at fixed parameter k). The problem with 
this family of analytical solutions is that, at large distances, all 
components (which have the same radial dependence for any /) 
decay too slowly: Br r~^,Be r~^,B(p r~^. Thus this con- 
figuration cannot be a solution for the whole space, as it would 
imply infinite magnetic energy in an infinite volume. Also, these 
solutions cannot be continuously matched with vacuum, because 
it would require that, at the same radius rout, ^cj^i^out) = and 
^riTout) ^ 0, a condition that cannot be satisfied because those 
two components are both proportional to the Yi{x) functions and 
have the same zeros. 



2.3. Self-similar models 



iLow & Loul (119901 followed by "Wolfson| (Il995h and other au- 
thors, studied a particular class of self-similar solutions to de- 
scribe the op ening of the solar coronal magnetic field lines due 
to their shear. [Thompson et"aD (I2QQ2 ) applied the same approach 
in the magnetar framework. Assuming a radial power-law form 
for the magnetic flux function and a radial dependence a oc 1/r, 



r = ro(y)Vou), 

I r-/ \\l/p I l-^l 

r \ro 



the enclosed current function is 



m = k 



where 



i+i/P 



4(p+l) 



(20) 
(21) 

(22) 
(23) 



is the current normalization, related to the amount of current, like 
the dimensionless parameter kss- Equation (fTOb becomes a non- 
linear second-order difl'erential equation for the angular function 



(l-fi^)F''(ji)^p(p^l)F(ji) 



> + l 



FmFiMf^P (24) 



where primes stand for derivatives of F(fi) with respect to yu. The 
magnetic field is given by 



Bo 
2 



\(p+2) 



F(m) 



(25) 



(p+2) 



p F(Li)\F(u)\ 



Self — similar family 
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Fig. 2. Curve Io(p) describing the family of self-similar twisted 
dipoles. 



The three components of the magnetic field have the same 
radial dependence, so we are seeking for self-similar solutions. 
The maximum line twist, which is the one suff'ered by a polar 
line (r ^ 0), is called global twist, and it can be taken as the pa- 
rameter that uniquely describes the /-pole family of self- similar 
solutions. Due to self- similarity, Eq. (O for F ^ becomes very 
straightforward: 



A(p : 



2k 

p 



-hi Jo 1-yU 



Up 
2 



(26) 



Equation (l24l) has to be numerically solved by imposing 
three physical requirements: F(l) = F(-l) = (on the axis, 
only the radial field component is allowed), and f (1) = -2 (to 
fix the normalization). Given a value of kss, one can find an infi- 
nite number of solutions characterized by an eigenvalue p. Each 
solution represents a diff'erent multipole (except the dipole for 
which there are two solutions). Fixing kss = (i-e. /o = 0, no 
current) implies an integer value of so that we recover the vac- 
uum multipolar solution with F(ji) = (1 - ji^)dPp(ji)/dji where 
Pp is the p^^ Legendre polynomial. For each multipole, a unique 
relation kss(p) (or Io(p) once the values of and are fixed) 
defines the family of solutions. In Fig. [2] we show the curve of 
parameters we obtained for the twisted dipole family. Our re sults 
agree with the dipolar solutions of Thompson et al. j2QQ2h and 
the higher multipole solutions of Pavan et al. ( 2009). We have 
used these models in the numerical code (section Sect. [3]) for 
testing purposes and for the sake of comparison with other nu- 
merical solutions. 

However, all these solutions are of limited generality because 
they have a defined symmetry with respect to the equator, and a 
linear combination of solutions for diff'erent multipoles is not a 
new solution, due to the nonlinear character of the problem. 

2.4. General axisymmetric solutions 

The righthand side depends on the functional form of a(T) and 
is responsible for the coupling between diff'erent multipoles. We 
have explored several possible analytical forms of a(T). Aiming 
at more general solutions, we explore several possible analytical 
forms of a(T) in Eq. ([T5]) . In search of semianalytical solutions, 
we consider power laws of the form 
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with k the dimensionless parameter regulating the amount of cur- 
rents, and Fo given by Eq. ([T2l) . 

Obviously the trivial choice = gives the vacuum (i.e. un- 
twisted) solutions. The choice of a constant has already been 
considered in Sect. 12.21 Some choices with ^ < 0, such as 
q = -1/2, present mathematical problems and are definitely not 
physical, and high values of q make the problem unsolvable with 
semi-analytical methods. Therefore, only a few simple choices 
may be useful because of their simplicity. An interesting case is 
q = 1/2. Assuming that F > in the domain of integration, this 
particular choice leads to 



X B.y is a normalization constant that sets the 



F 

V sin ^ \ Fo 



3/2 



(28) 



and the following ODE for each radial function // = a/(r)r/r^: 



/(/ + = - Z f^^'-^M^^Slmn. (29) 

^ ' m,n-l 



where each dimensionless Gaunt factor gimn involves the integral 
of the product of three Legendre polynomials (see Appendix |B] 
for a detailed derivation). Numerically it is ~ 0(1). Therefore, 
the source term couples each multipole / with an infinite number 
of multipoles m and n. One can solve the problem by truncating 
the series and allowing only for a finite number of multipoles. 
This approach is more general than other simple cases, but it has 
many more degrees of freedom. 

3. Numerical solutions 

In the previous section we have discussed some analytical and 
semi-analytical solutions that share the same drawback: the ar- 
bitrary choice of the enclosed current function /(F) or, equiv- 
alently, a(T). Some of these solutions are nonphysical, in the 
sense that they can neither be extended to infinity nor matched 
to vacuum solutions. All these limitations make the (semi- 
)analytical approach insufficient for general purposes, because 
we have no physical argument for preferring one particular form 
of the current over another. The alternative is to find numerical 
solutions of the nonlinear, force-free equations describing an NS 
magneto sphere. We have three main reasons for working in this 
direction. First, we expect these solutions to be more general and 
in some cases very diff'erent from the semi-analytical ones. By 
studying new numerical equilibrium solutions, we hope to gain 
insight into the form of the most realistic choices of enclosed 
current functions. Second, a reliable numerical code could in 
principle be extended, if adding the rotationally-induced elec- 
tric field, to build a consistent corotating magnetosphere model 
from which the distribution of charge density can be computed 
directly. Third, the ability to build magneto spheric solutions for 
any prescribed form of the magnetic field at the surface is very 
useful for future studies of the dynamics of the magnetosphere 
coupled with the physics of the crust. 

3.1. The magneto-frictional method 

For simplicity, we begin with the two-dimensional nonrotating 
case, assuming that the force-free approximation is valid up to 
an outer radius rout, which reduces the problem to finding so- 
lutions of E g. ([TOl) . In th e so-called magneto-frictional method 
(I Yan g et aH, 1986, Roumeliotis et al., 1994), one begins with an 
initially non-force-free configuration and defines a fictitious ve- 
locity field proportional to the Lorentz force, Vf = vj x B/B^, 



where / 

time unit in the induction equation. Hereafter we use v 
results in a fictitious electric field: 



-VfXB. 



l.This 



(30) 



This fictitious electric field enters into the induction equation as 
a frictional term that forces the solution to relax to a force-free 
configuration. 

In the original method, iRoumeliotis et al.l (1 19941) write the 
magnetic field as 



5 = VF X VyS , 



(31) 



where F is the magnetic flux function defined in Eq. ([3]), and 
j3 = (p- y(r^ 0)^ where y{r, 6) is a scalar potential related with the 
freedom in the choice of the form of the enclosed current. The 
induction equation becomes a system of two advection equations 
for those two functions: 



dtY -h V/ • VF = 



(32) 



A static solution is achieved if and only if v/ = 0, since the 
velocity field, VF, and are all orthogonal to the magnetic 
field by definition. 

We apply the same idea but, instead of evolving the functions 
F and we evolve the magnetic field components directly by 
solving the induction equation 

dtB = -VxEf, (33) 

where the fictitious electric field can be written as 

Ef = J-(JB)B/B\ (34) 

where ^/ is a measure of the deviation from the force-free con- 
dition, because J \\ B is accomplished if and only if = 0. The 
main reason for solving the induction equation instead of Eqs. 
(l32l) is to allow for future extensions of the code by considering a 
real, rotationally-induced electric field. The disadvantage is that 
we have to be more careful when setting boundary conditions for 
the electric field, because we could converge to stationary solu- 
tions characterized by V x = 0, which are not necessarily 
force-free. 



3.2. Linear analysis of the magneto-frictional method 

We now consider a background, uniform magnetic field and 
a small perturbation SB oc e^i^-r-ajt) ^ the linear regime, the 
equations read as 



6J = —ikxSB 
An 



dSB 



[(SJ X Bo) X Bo] 



(35) 



= -icoSB = -V X SEf . 



Explicitly, the last equation can be written as 

_ i^sB = ^^^[6B • (k X Bo)] - k^6B . 
c Bl 



(36) 



If the perturbed current is orthogonal to the background mag- 
netic field (either longitudinal perturbations 6B \\ Bq or trans- 
verse perturbations with k \\ Bq), the first term on the righthand 
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Fig. 3. Location of the variables in the numerical grid. 



side of Eq. ([36b vanishes and the dispersion relation is purely 
dissipative: 



oj = -i — k . 
An 



(37) 



Any perturbation of this type will be dissipated on a timescale 
oc k~^. In contrast, for transverse perturbations with both SB and 
k orthogonal to ^o, the current is parallel to the background 
magnetic field, and the two terms in the righthand side of Eq. 
([36b cancel out, so that the perturbation does not evolve (a neu- 
tral mode with 6l) = 0). Therefore, the magneto-frictional method 
is designed to dissipate all induced currents nonparallel to the 
magnetic field but allows for stationary solutions with currents 
parallel to the field. Since the largest lengthscale in our problem 
is set by the size of the numerical domain, Amax^^out^ the typical 
diff'usion timescale on which we expect to converge to a force- 
free solution is tdif oc r^^^. 

3.3. The numerical method 

We work in spherical coordinates (r, 6, (p) under the assumption 
of axial symmetry. We emplo y a fully explicit fin ite diff'erence 
time domain (FDTD) method (ITaflove & Brodwin[ [T975) with a 
numerical grid equally spaced in 9 and logarithmic in the radial 
direction, unless the outer radius is very close to the surface of 
the star (located at r = 1), in which case we employ a linearly 
spaced radial grid. Our typical resolution varies between 30-200 
points in the radial direction and 30-100 points in the angular 
direction. At each node (^/, r^), we define all components of 5^^'^^ 
and consider the three surfaces 5"^^'^^ centered on this node and 
normal to the unit vectors h = r,6, 0. The elementary areas are 
calculated by the following exact formulas: 

Sr = \ r^i sinOdcpde = 27rr^(cos6/,_i - cos^,+iP8) 
Jo J0,_i 

X27T r>rj+i 
j r sin Oi d0 dr = 7i{r]^^ - r]_^) sin Oi (39) 

r/+i rrj+i 
rAOdr = {Om - - r^j.,)l2 . (40) 



Since the the righthand side of Eq. ([33b contains V x it is 
useful to apply Stokes' theorem as follows. The magnetic flux 
across S^'^^ is approximated by 



and the magnetic fluxes are advanced in time using 

& Ef . d/, 



^t'^\t -h AO = ^t'^\t) - At 



(41) 



(42) 



where the numerical circulation of the electric field along the 
edges of the surface 5„ is approximated by using the values of 
Ef in the middle of the edges, whose lengths are 



(43) 
(44) 
(45) 



Explicitly, the circulation of £/ along the edge of each surface 



S^-^ is 



^5' 



E'dl = 



E-dl 



E'dl. 



I-J f if \ J-J f iy 



(46) 
(47) 

(48) 



In Fig. [3] we show the location of the variables needed for the 
time advance of B^]~^'^\ B 
We also make use of S 
components at each node: 



time advance of B^^ (^j^^d), and B^^'^^ (green). 

We also make use of Stokes' theorem to calculate the current 



Then, the values of /^^'^^ and B^^'^^ directly provide E^^'^\ as de- 
fined by Eq. ([34b . For each cell (/, j) the local divergence of B 
is defined as the net magnetic flux flowing across the surfaces 
divided by the cell volume (fluxes through the toroidal surfaces 
S(p cancel due to axial symmetry): 

(V . Bf^j^ = [O^^"'^"^'^ - ^^^-'^ + ^^^'^^^ - ^l'^^^]/V^'^j\ (49) 
with 

r>2n r^Oi+i r^rj+i 

y^'^j^ = r2sin6/d0d6/dr 

Jo JOi-i Jr^_i 



In 



= —(cosOi-i - cos6//+i)(rJ^i - rj.^) , 



(50) 



The numerical method ensures that the local divergence is con- 
served to machine accuracy by construction. As a matter of fact, 
Eqs. ([42b and ([49b imply 

d(V • Bf^j^ 



dt 





•d/-h 


d) E 


•dl 












d/-F ( 


£ E 


dl 











(51) 
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According to Eqs. (06]), the last equation is written as a sum 
of toroidal elements ^^^Z^. They are evaluated twice at each of 
the four surrounding edges (with center ^/±i,ry±i), with oppo- 
site sign, so they cancel exactly. Therefore, provided there is an 
initial divergence-less field, it is guaranteed that the numerical 
divergence defined by Eq. (|49l ) remains zero to round-off' level at 
each time step. 

We must mention that we also tried a method with a stag- 
gered grid (Yee, 1966) , in which each ^-component of the mag- 
netic field is defined only at the center of the normal surface, 
Sn'^\ while the electric field components are defined in the mid- 
dle of its delimiting edges. Methods based on staggered grid are 
well- suited to solving Maxwell's equations, because it provides 
a natural way to time-advance one field by means of the circula- 
tion of the other one. However, we are not dealing with the true 
Maxwell's equations, but rather with an artificial electric field, 
Eq. (l34l) . Evaluating the dot product J • B requires the interpo- 
lation between two or four values of several of the six mutually 
displaced components. Considering the red components in Fig. 
[3l for instance, the calculation of E^'^^'^^ also requires B^^'^^'^\ 
which would be not defined at that location. The unavoidable 
interpolation errors prevent the code from completely relax to 
Ef = 0, except in the trivial case of untwisted configurations. 
For this reason, we decided to work with a standard grid. 



3.4. Boundary conditions 

At the polar axis, we impose the vanishing of all angular compo- 
nents of magnetic field and currents: Eq = B^p = Jq = J^p = Eq = 
E(p = . At the surface, we have to fix the magnetic field com- 
ponents as provided by some interior solution. However, an ar- 
bitrary choice of poloidal and toroidal fields may not be compat- 
ible with a force-free solution. We decided to impose ^/ = 0, at 
the surface, which is equivalent to keeping the value of the radial 
component fixed at the surface, Br(l,0) and therefore to fixing 
the angular dependence of the magnetic flux function, r(l, ^). As 
a consequence, Bq and B(p are allowed to vary on the first radial 
grid point. 

The external boundary is set at r = rout- We have explored 
two diff'erent boundary conditions: Ef(r > rout) = and the 
continuous matching to external vacuum solutions. The first 
choice is equivalent to fix the radial component Br(rout, 0), while 
allowing for Bq and B^p to evolve. Coupling to vacuum solu- 
tions can be done using the spectral Legendre decomposition 
of the radial field at the outer surface, or more precisely, of the 
magnetic flux function, in terms of which the vacuum bound- 
ary condition for eac h multipole is easily imposed (e.g., as in 
iPon s & Geppert', 2007). Then we can reconstruct the angular de- 
pendence of Bo(rout, 0). The vacuum region is characterized by 
B(p = and the absence of currents or fictitious electric fields. 
This implies that a current sheet Je(rout) ^ is needed to ensure 
current conservation. 

If we choose Ef = as outer boundary condition, the 
code can actually converge to ^/ = at a round-off' level, be- 
cause mathematically this is the only solution compatible with 
V X Ef = 0. The price to pay is a forced matching of the inner 
solution with a fixed value of Br(rout)- In contrast, if we cou- 
ple with vacuum, there is no guarantee that the final solution is 
Ef = everywhere. We discuss below the inffuence of the dif- 
ferent boundary conditions on the results. 



^ out 




He 


tdis 


5 


30 


30 


41 


5 


50 


100 


35 


10 


30 


30 


183 


10 


50 


30 


162 


100 


50 


30 


180 X 102 



Table 1. Time needed to dissipate the numerical currents of the 
vacuum dipole. 



3.5. Convergence criterion and tests 

Since the magneto-frictional method is based on introducing an 
artificial, viscous electric field that drives an arbitrary initial con- 
figuration into a force-free state, we need a convergence criterion 
to decide when our solutions are acceptable. For that purpose, we 
keep track during the run of the following quantities: 

- Volume-integrated magnetic energy (total and contribution 
from the toroidal field) 

r B^ rBl 

St = J ydV, 8t^ = j -fdV. 

- Volume-integrated energy stored in the fictitious electric 
field 

- Total volume-integrated helicity 

^= f B^A^dV 
Jv 

(see Appendix lAl for a discussion about this definition). 

- Volume-integrated absolute value of the divergence of both, 
B and /. These two quantities are expected to vanish at 
round-off' level by construction. 

- An average throughout the entire magnetosphere volume of 
the local angle between current and magnetic fielc^ 

. 2 - Z/'sin^y/ J 

sm,= ^^ = ^^. (52) 

where the sum is performed over each node (/, j). 

- A last important cross-check is the consistency of the func- 
tions /(F) and a(T). First, we check that, for each n- 
component, the three functions Qf„(r, 0) = 47iJn/cBn(r, 6) are 
the same. Second, at each radius or angle, the relation Q has 
to be satisfied. 

Hereafter, we show the electromagnetic quantities in units 
of ^★^ so the magnetic field B scales with ^o, the current 
density / with cBo/ri,, the enclosed current / with cBor^,, and 
the magnetic ffux function F with Tq = Borl/2.To test our code 
and to fix our convergence criteria for the realistic models, we 
performed a battery of tests. In the first basic test we considered 
the analytical vacuum dipole with = Bosm6/2r^, B(p = 
and checked the ability of the code to maintain this solution. 
Due to the discretization errors, a little numerical toroidal cur- 
rent appears, and consequently a nonvanishing toroidal fictitious 
electric field. These errors are small (A8b/Sb < 1%), but it is 

^ This average weighted with avoids numerical problems in re- 
gions where the numerical value of the current is very low and the angle 
sin^ ri - E ■ J I is numerically ill-defined. 
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Fig. 4. Influence of the outer boundary condition and the location of the external boundary on the evolution of the ratio Bel^b(p (left) 
and 7) (right). 



Current density 



Current density 



Current density 



Current density 



12 3 



Fig. 5. Current density distribution / near the surface for solutions obtained with the same initial data and boundary conditions, 
Ef = 0, but varying rout = 5, 10, 50, 100 (left to right). 



interesting to see how long it takes to dissipate the perturbation 
to obtain 6^ = to machine accuracy. In Table [T] we show the 
results with difl'erent resolutions and values of outer radius, al- 
ways keeping the timestep close to the maximum value allowed 
by the Courant condition. The expected behavior tdis ~ r^^t 
obtained, with the constant of proportionality depending on the 
grid resolution, which afl'ects the strength of the initial numerical 
current. 

The second, less trivial test is provided by the analytical 
models discussed above: the Bessel and self-similar solutions 
(Sect. 12.21 Sect. 12.31) . Again we began with an initial model con- 
sisting of a known solution, and let the system dissipate the cur- 
rents that come from discretization errors. To obtain the initial 
models, the self- similar solutions require the numerical resolu- 
tion of the nonlinear ODE (1241) , while the analytical Bessel solu- 
tions are directly implemented. We tried difl'erent parameters for 
the Bessel solutions (varying k) and self-similar models (vary- 
ing the multipole index and the global twist). For every model 
tested with analytical solutions, we observe a very slight nu- 
merical readjustment of the configuration and the code rapidly 
reaches the relaxed state, with relative changes in St, Sbcp^ ^ 
less than ~ 1%. We had to pay special attention to work with 
suflftcient radial resolution in the case of highly twisted Bessel 
models ^ > 1, due to their oscillatory radial dependence. For low 
resolution, the code may find, after a large scale reconfiguration, 
a completely diff'erent solution with smaller which is numeri- 



cally more stable. If the resolution is high enough, all analytical 
solutions are found to be stable. 

We have also studied the evolution of vacuum dipolar solu- 
tions with an additional toroidal field for diff'erent values of Vout- 
In this case the initial currents are due not only to discretiza- 
tion errors, but also to an inconsistency in the initial model. 
Moreover, the mean angle defined in Eq. (l52l) initially has a fi- 
nite value 7) in. In Fig. |4] we show how some convergence mon- 
itors evolve as a function of time for three cases with an initial 
toroidal field of the form = OABosinO/P {fUn = 15.1°), but 
diff'erent external boundary conditions: matching with vacuum 
at rout = 10 or imposing = at rout = 10 or 100. For compar- 
ison, we also show results for two analytical solutions: a Bessel 
solution and a twisted self- similar model with the same helicity. 

Finally, we tested a vacuum dipole perturbed by a weak 
toroidal field. This is a case of physical interest for quasi pe- 
riodic oscill ations of magnetars , as discussed in Ti mokhin et al.l 
(l2008h and iGabler et all(l201ll) . Given a background poloidal 
field described by F, we chose an arbitrary functional form /(F) 
and built the toroidal field according to Eq. (O. As expected, 
the perturbed configuration is stable and the stationary solu- 
tion is rapidly reached after a small readjustment. Typically, for 
msix(B(p) = 0.1 we have fjin ~ 0(1°) and changes AB/B ~ 
1%. 

In general, the magnetic energy is not conserved, since the 
system has to dissipate part of the current to reach a force-free 
configuration. This eff'ect is more evident for initial configu- 
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0.0 0.2 0.4 0.6 o.a 

Fig. 6. Enclosed current I(T) for different boundary conditions 
= or coupled with vacuum) at rout = 5, 10, 50, or 100. 



Model 


kfor 


ginir) 




lin [deg] 


A 


0.010 


r-' 


mi 6 


13.8 


B 


0.100 


r-3 


sin^ 


15.1 


C 


0.500 


r-3 


sin^ 


29.9 


D 


0.115 


r-3 


sin^6' 


10.7 


E 


0.230 


r-5 


sin^6' 


25.1 


F 


0.234 


^-3^-[(r-l)/0.5]2 


sin^6' 


23.3 


G 


0.113 


r-3 


sin^6' 


34.2 



Table 2. Parameters defining the initial toroidal field, as indi- 
cated in Eq. ([53]) . 



rations with high helicity. When the outer boundary condition 
Ef = is imposed, the helicity is conserved within a few per- 
cent, as expected (see Appendix lAl for the helicity conservation 
theorem), and both electric field and mean angle eventually van- 
ish (to machine accuracy). However, when rout is large, or when 
vacuum boundary conditions are imposed, configurations with 
high initial helicity take a much longer time to relax (see Fig. 
131). In all cases, the relaxation process is faster near the surface, 
where the configuration of the magnetosphere is more important 
for our purposes. 

On the basis of these results, our convergence criteria for ac- 
cepting that a configuration has reached a force-free state are 
hereafter 6^/5^0 < 10"^ and fj < 10"^ °, plus the require- 
ment that both quantities are monotonically decreasing with 
time. Some short, initial relaxation phase, in which some large- 
scale reconfiguration occurs is possible. We chose to compare 
the electric energy to the magnetic energy contribution from the 
toroidal field, which is much more restrictive than simply the 
ratio of electric to magnetic energy, especially for low helicity. 
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9{ 






Jr^ 


n 

F 








\rBnlr^^ 

LCUQ/ ' ★J 






SI 


0.21 


0.5 


1.8 X 10""^ 


0.061 


0.97 


S2 


1.11 


1.6 


8.1 X 10~^ 


0.15 


0.69 




0.021 


0.12 


1.2 X 10~^ 


no'so 

\J .\J\J U\J 


1.45 


D 


91 


1 9 






1 40 


c 


1.11 


6.5 


4.8 X 10-2 






D 


0.21 


0.7 


1.5 X 10-2 


0.056 


1.14 


E 


0.21 


0.4 


4.2 X 10-2 


0.106 


0.50 


F 


0.21 


0.4 


4.3 X 10-2 






G 


0.21 


0.3 


2.0 X 10-2 


0.038 


2.22 



Table 3. Comparison between two self-similar solutions (S1-S2) 
and our numerical solutions (A-G). 



4. Results 

With the numerical code described above, we can obtain very 
general solutions of force-free, twisted magnetospheres. We dis- 
cuss separately the influence of the following relevant parame- 
ters: 

- the location of the outer radius rout in combination with the 
external boundary condition; 

- angular and radial dependence of the initial toroidal field; 

- initial twist and helicity, fixed by the functional form and the 
strength of the initial toroidal field; 

- the geometry of the initial poloidal field. 

In Fig. [5] we compare the distribution of currents in a solution 
obtained by imposing = at rout = 5, 10,50, or 100. In 
all cases the initial poloidal component is a vacuum dipole so- 
lution with A(p = BosinO/r^ and a toroidal field of the form 
B(p = OABq sin 6/r^. We observe that, near the axis, the solutions 
are clearly aff'ected by the location of the external boundary, if 
it is not far enough from the surface (rout ^ 10). In such a case, 
the influence of the external boundary is important, and it intro- 
duces artificial features, although the current distribution in the 
equatorial region is similar in all cases. The final configurations 
become almost indistinguishable when rout = 50 or 100. Taking 
rout ^ 100, we ensure that the numerical noise caused by the in- 
teraction with the external boundary is negligible. By neglecting 
the contribution from the open field lines, for which the twist is 
ill-defined, the global amount of twist is similar in all models 
(~ 1.2 rad). 

The function /(F) for the same four cases is shown in Fig. 
[6l together with two more cases corresponding to rout = 5 and 
rout = 100 but replacing the external boundary condition Ef = 
by a smooth matching with vacuum solutions. The use of a dif- 
ferent boundary condition afl'ects the final solution only for low 
values of rout- Matching with external vacuum implies that no 
currents can flow through the boundary, so that /(F) = a(T) = 
along every field line crossing the outer boundary. As a conse- 
quence, a force-free configuration, coupled with a vacuum, will 
be characterized by a plateaux /(F) = for F < F^, with F^ label- 
ing the first closed field line. This means that a bundle of open 
field lines is potential. The fraction of open field lines (the length 
of the plateaux for low F) is large only for low values of rout- 
For F > 0.2 (equatorial region), all curves coincide. Increasing 
rout > 100 further has no visible eff'ect on the models with this 
helicity. For models with higher helicity, the interaction with the 
boundary becomes more important, and rout needs to be accord- 
ingly increased to minimize the boundary eff'ects. 
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Fig. 7. Self-similar model S2. From left to right: poloidal magnetic field lines (white) and strength of the toroidal component (colored 
logarithmic scale); current density distribution / in the near region (r < 5) (colored linear scale); angular profiles of the magnetic 
field components close to the star surface, and angular profiles of the current density. 
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Fig. 8. Same as Fig.|7]for model A. 
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Fig. 9. Same as Fig. [7] for model B. 



The next step is to explore the influence of the strength and 
form of the initial toroidal field. Table [2] summarizes the param- 
eters of the initial models employed and the initial mean angle 
between B and /, Eq. (l52l) . Models A-F are obtained starting 
from a dipolar poloidal component with strength at pole and 
a toroidal component given by the general form 



The initial configuration of model G is asymmetric with respect 
to the equator: it is a superposition of dipolar and quadrupolar 
poloidal components plus a toroidal component, as follows: 



^0 sin 6 Bq sin 6 cos 6 
— ~ I" ~i ^; 



(55) 



B'^ = korBo gin(r)Me) . 



(53) 



B^ = OAB^ 



sin 6 



0- 



The angular part is chosen to be of the form fin(0) = sin^ 0, with 
d being a positive integer. The radial dependence of models A to 
E is a power law, gin(r) = r~\ while in model F we use a rapidly 
decaying function: 



- ^-3^-[(r-l)/0.5]2 



n(r) = r e 



(54) 



We fix rout = 100 and the external boundary condition to 
Ef(rout) = for all models, except for model C (highest he- 
licity), where rout = 500 to reduce the influence of boundary, as 
discussed before. 

Table [3] shows the features of our final configurations: he- 
licity, maximum line twist, maximum value of current density. 
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Fig. 10. Same as Fig.|7]for model C. 
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Fig. 11. Same as Fig. [7] for model D. 
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Fig. 12. Same as Fig.|7]for model E. 
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Fig. 13. Same as Fig.|7]for model F. 
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Fig. 14. Same as Fig.|7]for model G. 

and parameters of /(F) = /o(F/Fo)^^^^^. We also include (model 
SI and S2) two self-similar solutions (Sect. [231) of similar he- 
licity, where all components have the same radial dependence 
Bi ~ y-{p+^) ^ jj^ )^{^ case current function is analytical, 
/ = /o(F/Fo)^^^^^ with (/?, /o) belonging to the family of solu- 
tions in Fig.[2j In the other cases /q and p are obtained with fits 
to the numerical function. 

The final geometry of the magnetic field and currents for all 
these models is shown in Figs. [7] to [141 For ktor ^0.1 (mod- 
els A and B), the initial poloidal field remains almost unaltered, 
and the behavior of the solutions is nearly linear. Toroidal field 
strength, helicity, current density /, enclosed current function 
/(F), and global twist scale linearly with ktor- In contrast, for 
ktor ^0.1, the high initial helicity results in a larger twist angle, 
up to several radians, which in turn corresponds to a highly de- 
formed poloidal field. A direct comparison of models A, B and 
C, which diff'er from each other only in the strength of the initial 
toroidal component, illustrates this eff'ect: models A and B have 
the same shape with just a diff'erent scale factor, but model C is 
qualitatively diff'erent. 

The general features in the low ktor models do not diff'er 
much from self- similar solutions, because the initial conditions 
were close to a slightly twisted dipole. In self-similar models, 
two features are the absence of radial currents on the axis and 
a higher concentration of currents around the equatorial plane. 
Conversely, in our models, currents are spread more over the an- 
gular direction, and we allow for the existence of radial current 
on the axis. As a consequence, comparing numerical solutions 
A and B with self- similar models with comparable helicity, the 
former reaches lower maximum values of current density with a 
higher global twist. We also note that in the most extreme case 
(model C) the angular dependence of the toroidal field and ra- 
dial currents becomes steeper0 It is also interesting to compare 
models C and S2 (Figs. [TOl and [71). Both have a similar helicity, 
but the global twist is higher in model C, while the maximum 
current density is higher in model S2. 

Comparing models B and D, which only diff'er in the angular 
dependence of the initial data and the normalization (fixed to ob- 
tain the same helicity), we find that the converged solutions are 
very similar, except near the axis where model D has no radial 
currents. The eff'ect of varying the initial radial dependence can 
be estimated by comparing model D to model E. It seems that 
in this case the final solution keeps memory of the initial model: 



^ In some cases it approaches the formation of a current sheet near the 
equatorial plane, and this introduces numerical noise that does not allow 
reaching a smooth solution or calculating the twist angle accurately. 
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the converged solution shows a toroidal field that decreases with 
distance faster in the model E than in model D. 

In model E a tiny toroidal field (note the color scale in the 
figures) appears near the axis. This is likely a numerical artifact 
that can be partly ascribed to the (narrow) bundle of lines that 
depart from polar region and interact with the outer boundary. 
As a matter of fact, these structures are stronger for low values 
of Tout, as already shown in Fig. [3 Moreover, the numerical dis- 
sipation of the current is slower near the axis and longer runs are 
needed to reach more restrictive convergence criteria. 

The diff'erent radial dependences of the magnetic field com- 
ponents in models B, E and F, together with SI, are shown in 
Fig. [151 The self- similar solution has the same radial dependence 
y-{p+2) ^Qj. tjij-ee components, but in the numerical solutions 
the toroidal field can decrease faster (model E and F) or slower 
(model B) than the poloidal components. Furthermore, the radial 
behavior may depend on the magnetic colatitude 6, too. In addi- 
tion, the radial dependence of the poloidal components is close 
to the power law r"^ when the twist is low, but it may signifi- 
cantly deviate from a power law for the more complex models. 
The departure from self- similar solutions is likely to have a vis- 
ible eff'ect on the observed spectrum, as discussed in more detail 
in the next section. 

A diff'erent solution that can only be found numerically is 
model G, the asymmetric configuration. Most of the current den- 
sity is concentrated at a high latitude 6m. 

An interesting case (not shown) consists of an initial model 
with toroidal and poloidal fields of opposite parity (e.g. a dipolar 
poloidal field plus a quadrupolar toroidal field). In this case, the 
total initial helicity is zero and, since this is a conserved quantity, 
the current is dissipated and a potential solution is found by the 
numerical code. This is consistent with the fact that the numer- 
ical evolution converges towards the most trivial solution with 
the same helicity. 

Finally, to better understand the diff'erences among models, 
a comparison between the enclosed current function /(F) is very 
helpful, as shown in Fig. [161 In principle, if we know this func- 
tion or an approximation fits the results of numerical simula- 
tions, the pulsar equation (Eq. [71) can be solved and the magne- 
tospheric structure can be determined. In our models A, B, D, 
E, G, the current /(F) is monotonic and can be well fitted by a 
single power law /o(F/Fo)^^^^^, but the values obtained for the 
parameter p are not consistent with the value of /q describing the 
self- similar family of solutions (Fig. [2l), as expected. As a mat- 
ter of fact, for models A, B, D, the values of p are greater than 
1 (the self- similar dipolar family is described by p e [0, 1]). In 
contrast, model E lies quite close to the self-similar solution. The 
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Fig. 15. Radial profiles of the magnetic field components for 
models SI, B (top), E and F (bottom). We show Br on the axis 
and Bo, B^p at the equator for r g [1, 30]. 



enclosed current of models C and E varies more rapidly, describ- 
ing the concentration of currents in a smaller angular region. For 
some models, a power-law fit to /(F) simply does not work. In 
model G, even if the symmetry with respect to the equatorial 
plane is broken, the resulting enclosed current /(F) is not very 
diff'erent from models A, B, D, and SI (after rescaling accord- 
ingly to the factor ktor)- The diff'erence is the mapping between 
F and the surface footprints. In all symmetric models, the maxi- 
mum magnetic flux (proportional to F) is located at the equator, 
while it corresponds to ^ n/2 for model G. Depending on the 
viewing angle, this has a strong imprint on the processed spectra. 

4.1. Resonant optical depth 

A prime candidate mechanism for generating the hard tail com- 
ponent observed in magnetar X-ray spectra is the resonant 
Compton upscattering of photons in an external magnetic field. 
We denote the cyclotron frequency by ojb = ZeB/mc, where Ze 
and m are the charge and mass of the particle. For electrons, this 
corresponds to an energy Hcob = 11.6(5/10^^ G) keV. In the X- 
ray band, the plasma frequency is much lower than the photon 
frequency, thus taking the polarization vector orthogonal to the 
propagat ion direction k (sem itransverse approximation) works 
weU (see lCanuto et al.1 (|l97lh and Ventura (1979)). The normal 
modes of propagation are in general elliptically polarized and the 
degree of ellipticity depends on the ratio oj/ojb and 6kB, defined 
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Fig. 16. /(F) for diff'erent models. The curve of model A is mag- 
nified by a factor of 10. 



as the the angle between k and 5. It is common to introduce 
the ordinary (O) and extraordinary (X) modes. In the limit of 
propagation perpendicular to the magnetic field (6kB = 7r/2), the 
modes are linearly polarized, with the polarization vector paral- 
lel (0-mode) or orthogonal (X-mode) to B. For parallel propa- 
gation (6kB = 0), we recover the circularly polarized modes. 

The cross section of O and E modes depends in a nontrivial 
way on 6kB, and on the ratio oj/ojb. The cross section of the O- 
mode is close to the Thomson cross section ctt for OkB near to 
7r/2, it scales with sin^ 6kB if OJ <^ ojb, and it does not depend 
dramatically on the frequency. Conversely, for low photon ener- 
gies fico <^fi(jjB, the cross section of the X-mode is strongly sup- 
pressed due to the reduced mobility of charged particles across 
magnetic field lines. Moreover, at photo n frequency oj = ojb, 
the X-mode becomes resonant (IVenturaL 1979). If we neglect 
the thermal velocity of the particles (cold plasma approximation) 
and the natural width of resonance due to the radiation-damping 
force, the resonant cross section can be approximated by a delta 
function: 



2 2 (^^) 
(Tres =71 (I -\- COS OkB) ^(OJ - COb) 



(56) 



where the factor (1 -h cos^ OkB) arises from the assumption of 
unpolarized light. The first attempts to co nsider this process 
quantitatively in magnetars are presented by lLyutikov & Gavriill 
(2006). They study a simplified, semi-analytical 1-D model by 
assuming that seed photons are radially emitted from the NS sur- 
face with a blackbody spectrum (BB) and the resonant Compton 
scattering occurs in a thin, plane parallel magneto spheric slab. 
By assuming a bulk velocity of electrons, but neglecting all ef- 
fects of their recoil (Thomson limit Hcob <^ rUeC^/y, where y 
is the Lorentz factor of electrons), it can be est imated how a 
BB spectra is aff'ected. iLyutikov & Gavriij (I2006h find an aver- 
age upscattering for the transmitted radiation (forward scattered 
photons), while the mean energy of the reflected radiation re- 
mains the same. Later, iNobili et al.l (l2008attbl) predicted spectra 
obtained via a new 3-D Monte Carlo code, accounting for polar- 
ization, QED efl'ects, and relativistic cross sections. Although in- 
cluding our magnetospheric models in such sophisticated Monte 
Carlo simulations is beyond the scope of this paper, we can es- 
timate some of the possible efl'ects by looking at the behavior of 
relevant quantities. 
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In particular, we can estimate the resonant optical depth, 
given by the integration of Eq. ([56b along a given line of sight of 
constant 6 (.Lvutikov & Gavriil . 2006): 



Tres(0) = 71 nzZe{\ + COS OkB) 



dB 
dr 



(57) 



where nz is the density of scatterer particles, and all the quanti- 
ties are evaluated at the resonant radius Vres- As the latter depends 
on the photon energy, the energy dependence of the optical depth 
is given basically by the local ratio (1 +cos^ 6kB)nzl\dBldr\ (pro- 
vided that Vres hes above the star surface). If we assume charge 
separation, the particle density is proportional to the current den- 
sity, / = KcnzZe. The dimensionless factor k depends locally on 
the plasma composition and on t he (bulk and/or thermal ) ve- 
locity of the particl es. Following [Thompson et al.l (l2002h and 
iNobili et al.1 (l2QQ8ah . we can (very roughly) estimate the opti- 
cal depth along a line of sight (assuming k constant along it), 
and considering only radially directed photons: 



dB_ 
~d^ 



(58) 



If we consider self-similar magnetosphere models under 
these naive assumptions, the optical depth becomes independent 
of Vres, hence, independent of where the scattering happens. This 
is because the local ratio (1 + cos^ OkB)J/\dB/dr\ is the same 
for each angle, since for every component /, Bi oc r"^^^, and 
Ji oc r"^^^ . Furthermore, the optical depth does not depend on 
the normalization of the magnetic field, ^o- 

In contrast, the optical depth in our numerical models de- 
pends on the photon energy, because the components of J and 
B cannot be described by the same power law. In Fig. [T7] we 
show the estimated resonant optical depth, Eq. ([58b , for our mod- 
els, compared with the self-similar model SI. We plot KTfgs as a 
function of the polar angle 0, for diff'erent energies of the seed 
photons, hoj = 0.5, 1, 2, 4, 8 keV, taking = 10^^ G. If the mag- 
netic field is predominantly poloidal (low helicity), the optical 
depth is roughly given by the ratio between toroidal field (which 
provides poloidal current) and poloidal field at the resonant ra- 
dius. Increasing the photon energy, the resonant radius will be 
closer to the surface, modifying the estimate of / and dB/dr 
which are involved in Eq. ([58b . Thus, if the toroidal field decays 
slower than the poloidal component, a higher photon energy im- 
plies a lower ratio J/\dB/dr\ at the resonant radius. Looking at 
the radial profiles (Fig. O, we can understand why the optical 
depth increases with the photon energy for models E and F, while 
for other models the optical depth is greater for soft photons. A 
more precise prediction of how the spectrum is modified when 
using one or another magnetosphere model requires more elab- 
orated calculations than those presented here. We also point out 
that, due to the linear relation between B and ojb, increasing the 
energy is equivalent to decreasing the normalization by the 
same factor. 



5. Conclusions 

The study of force-free magnetosphere models from an ana- 
lytical or semi-analytical point of view have led to remarkable 
results in a few cases. In the literature, the twisted configura- 
tions usually implemented in the modeling of synthe tic spectra 
are oft en based on the self- similar solutions by Tho mpson et al.l 
(I2OOI . which are re stricted to a twist ed dipole or to a sin- 
gle higher multipole dPavan et al.L 120091) . Since semi-analytical 



studies of a general combination of multipoles are more diffi- 
cult due to the nonlinear character of the equations, we faced 
the problem by constructing and thoroughly testing a numerical 
code that relaxes an arbitrary initial model to a force-free solu- 
tion with given boundary conditions. With our numerical simu- 
lations we can find general solutions of twisted magnetospheres 
with complex geometries. Our only input is the radial field at 
the surface, which can be completely general (i.e. provided by 
numerical simulations of the internal evolution). 

Our main conclusion is that the possible family of relaxed 
co nfigurations depends on th e initial data. As already pointed out 
by Contopoul os et aP ([201 1 ) in a different context, the solution 
is not unique, given the surface radial magnetic field and the 
outer boundary conditions. This reffects the freedom we have in 
choosing the enclosed current function /(F), for a given r(^) at 
the surface. Starting with two diff'erent configurations with the 
same value of volume-integrated helicity does not necessarily 
lead to the same final configuration. Thus we were able to find 
a variety of solutions, which can be qualitatively diff'erent from 
the self- similar models. 

Some of the new magnetosphere models that we found are 
characterized by a high degree of nonlinearity, twists up to a few 
radians, and a nontrivial functional form of the threading current 
/(F). These models establish a new basis for generalizing the 
study of radiation transfer in neutron star magnetosphe res, as al- 
ready extensively done with th e self-similar models (R ea et all 
I2OO8I iFernandez & Thompson! 120071 IZane et al.l |2009). In par- 
ticular, our approach provides a more general framework for 
studying how the geometry of the magnetic field and currents 
affects the output X-ray spectrum of a magnetar. If resonant 
Compton scattering is a crucial ingredient in forming the spectra, 
relaxing the self- similarity constraint in the models may have a 
significant eff'ect, as our preliminary estimates of optical depths 
have shown. This opens the possibility of implementing alterna- 
tive models in simulations of radiative transfer. In the future we 
aim at providing spectra to be compared with observations. 
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Appendix A: Magnetic helicity 

For a given magnetic field configuration, the usual definition of 
magnetic helicity 



(A • B)dV 



(A.1) 



is not unique due to the gauge freedom in the vector potential, 
A — > A + V*?, where is a scalar function. However, in axial 



symmetry (no (/"-dependence), the following definition of helic- 
ity 



AV 



is gauge-independent because 'H ^ H + J B^V^^ = 'H. We 
have employed this definition of helicity in numerical code. It 
can also be shown that this quantity is conserved during the evo- 
lution. Taking the time derivative (denoted by dt) and using the 
induction equation, we have 



= j A^{d,B^)dV + j (d,A^)B^dV = 
= - J A^^ • (V X E)dV - j E^B^dV = 
= j [V ■(A^^xE)-E-Vx(A^^dV 
= S {A^4> xE)-dS- f{E- B)dV . 

JdV J 



-f 



EMV = 



(A3) 



If the poloidal electric field vanishes at the boundaries, the total 
helicity is conserved. 

With the definition (IA.2b . the explicit expression of the mag- 
netic helicity for self-similar models (Sect. 12.3b is 

(A.4) 



n , k fw:^. 



16 



D. Vigano, J.A. Pons & J.A. Miralles: Force-free twisted magnetospheres of Neutron Stars 



and for a dipolar Bessel model (Sect. 12.21) we have 



7i(jcrsin^6/dy = -kB: 



These integrals can be evaluated analytically using Gaunfs for- 
mula ^^^^3^ 



f&x. 



(A.5) 



/3 



Appendix B: A nontrivial semi-analytical 
configuration 

Here we show the derivation of Eq. ([29]) in Sect. 12.41 Substituting 
a = -^ir/Fol^/^ in Eq. ([B]), with T given by Eq. leads to 



l(Ul)\d\rai(r)) ^ ^^ai(r) 



2(-y 



{s-b-w) _ 



X 



2/+ 1 



t=p ^ 

a-\- b -\- c 



(b + v)l(c + w)l(2s-2c)lsl 
(b - v)l(s - a)l(s - b)l(s - c)l(2s + 1)! 
(a-\-u-\- t)l(b -\-c -u-t)l 
■ u - t)\(b - c + w + t)\(c -w-t)l 



(B.7) 



3 J-i djd \To) 
The integral at the righthand side can be written as 



(B.l) 



I 



1 dPf 

[X 



1 



^r^sgnCDdyu = :^[Thgn(T)Pi(ji)]\ + 
d/i Yl 



p = max[0, c - b - u] 

q = min[b -\- c - u, a - u, c - w] . 

The integral J3 is non-zero if and only if both the following 
relations are satisfied: 

- b-c<a<b-\-c,8i triangular condition; 

- a-\- b -\- c is even. 



^ 



r— sgn(r)P/Ou)dyu + 

1 d/i 



dr. 



I d/j. 



For each fixed / in Eq. (IB. 61) , the conditions above are satisfied 
for an infinite number of pairs (m, n), with gimn having numerical 
values of ~ 0(1). Thus Eq. (IB.4b couples each /-multipole with 
where the third term comes from sgn'(jc) = 2S(x). The first and an infinite number of other multipoles. 



third terms in this equation are zero. Next, we assume by sim- 



It should be noted that this is valid only if F > in its whole 



pHcity that F > in the whole domain, so that sgn(F) = 1, and domain (r,yu) e [1, rm ax] x [ -1, 1 ], otherwise we cannot use the 

we define the dimensionless radial functions /^(r) = am(r)r/r^. formulas above (Eqs. (IB.5I) , (IB.7I) ). A non positive F would make 

Hereafter we drop the dependences on yu, r, and the integration the calculation of these factors much harder, 
limits for conciseness. Using the Legendre equation we obtain 
from Eq. (fTTT) 



dr 



— = -Fo^^(n+l)/,P. 



(B.2) 



and we can express the integral J in a more compact form 

I = 2Y^n{n + \)fnY,fm {l-fi^)Pi-^Pn= (B.3) 



m-l 



m,n-l 



2m(m + 1) 
2m + 1 



PjP. 



where we have also used the recurrence relations between 
Legendre polynomials. Finally, the ODE ([T5]) for each //(r) is 

= _ _ ^/„/„g,„„, (B.4) 

^ ' m,n=l 



where 

2 2/+ 1 m(m+ 1) ^ 

3 2m + 1 /(/ + 1) 



J PlP^-lPn- J PlPm^lPn 



The integral of the product of three Legendre polynomials can 
be expressed by Wigner 3 - 7 symbols: 



J PaPbPc = I 



a b c 




(B.5) 



Alternatively, we can express the factors above in terms of asso- 
ciated Legendre polynomials: 



glm 



2n{n+ 1) 
3 /(/+ 1) 



f p^- f pi pi pO 



(B.6) 



